# code for "Racial Hierarchy and Jurisdiction in U.S. Status of Forces Agreements"
# Freeman
# Security Studies 

setwd()        # set working directory

require(maps)
require(viridis)
library(tidyr)
library(plyr)
library(dplyr)
library(moments)
library(xtable)
library("dplyr")
library(ggplot2)
library(arm)
library(lmtest)
library(sandwich)

# Figure 1: World Map of U.S. SOFA Jurisdiction

SOFA_map_data <- read.csv("SOFA_map_data.csv")  # load data for world map of SOFA jurisdiction from the CSV file into DataFrame. 

theme_set(
  theme_void()
)
map.juris.data <- SOFA_map_data        
map.juris.data <- map.juris.data %>%
  dplyr::select(StateNme, jurisdiction) %>%                      # select the two columns of interest: SOFA jurisdiction and State name
  dplyr::rename(region = StateNme, value = jurisdiction)   # rename columns

map.juris.data$region <- as.character(map.juris.data$region)                           
world_map = map_data("world")                         # retrieve the world map data
map.juris.map <- full_join(map.juris.data, world_map, by = "region")

ggplot(map.juris.map, aes(x=long, y=lat, group = group, fill = value)) +
  geom_polygon(colour = "white") +
  scale_fill_continuous(low = "orange",
                        high = "light blue",
                        guide= "colorbar")       # produces gradient scale legend. Instead, I construct legend separately within the .jiff file to represent geographic data on map, i.e. jurisdiction types (see figure in paper).   

# Figure 2: Interactions between racial hierarchy, extraterritoriality, and SOFA jurisdiction

# I construct this framework in Microsoft PowerPoint

# Figure 3: U.S. SOFAs by Jurisdiction, 1951-2020

# I construct this figure in Microsoft PowerPoint, referring to data: 

SOFA_overtime <- read.csv("SOFA_jurisdiction_overtime.csv")  # load data for SOFA jurisdiction over time from the CSV file into DataFrame. 
agg <- aggregate(cbind(sofa) ~ jurisdiction + year, data=SOFA_overtime, FUN = sum,na.rm=TRUE) # SOFA grouped by jurisdiction and year 


######################################## Models ############################################

SOFA_data <- read.csv("SOFA_data_2023.csv")

# Table 1: Descriptive Statistics

descriptive_stats <- SOFA_data %>% 
  dplyr::select(jurisdiction, nonwhite, cold_war, sharia, nato,
                veto, rule_law, judicial_ind, capability, affinity, troops) 

descriptive_stats <- na.omit(descriptive_stats)

# Prepare data 

SOFA_data [] <- lapply(SOFA_data, function(x) {
  if(is.factor(x)) as.numeric(as.character(x)) else x
})

skewness(SOFA_data$capability, na.rm = TRUE)            # fix high positive skew for capability indicator
SOFA_data$capability <- log10(SOFA_data$capability)

SOFA_data <- as.data.frame(SOFA_data)

# select variables for model

SOFA_models <- SOFA_data %>% dplyr::select(ccode, year, jurisdiction, nonwhite, cold_war, sharia, nato,
                                           veto, rule_law, judicial_ind, capability, affinity, troops)


SOFA_models <- na.omit(SOFA_models)

count(SOFA_models, 'jurisdiction') # exclusive = 40; concurrent = 143

SOFA_models$year <- as.factor(SOFA_models$year)
sapply(SOFA_models, class)

# Table 2: Non-White Majority Host on U.S. SOFA Jurisdiction

# race regression (rare logit) - column 1

race_rare <- bayesglm(jurisdiction ~ nonwhite, 
                      data = SOFA_models, family=binomial)
summary(race_rare)

coef_race_rare <- (coeftest(race_rare, vcov=function(x)
  vcovHC(x, method= "arellano", cluster="group", type="HC1")))

# law regression (rare logit) - column 2

law_rare <- bayesglm(jurisdiction ~ rule_law, 
                     data = SOFA_models, family=binomial)
summary(law_rare)

coef_law_rare <- (coeftest(law_rare, vcov=function(x)
  vcovHC(x, method= "arellano", cluster="group", type="HC1")))

# cold war regression (rare logit) - column 3

cold_rare <- bayesglm(jurisdiction ~ cold_war, 
                      data = SOFA_models, family=binomial)
summary(cold_rare)

coef_cold_rare <- (coeftest(cold_rare, vcov=function(x)
  vcovHC(x, method= "arellano", cluster="group", type="HC1")))

# all variables regression (rare logit) - column 4

main_rare <- bayesglm(jurisdiction ~ nonwhite + rule_law +
                            sharia + capability + troops + 
                            nato + affinity + veto + cold_war, 
                          data = SOFA_models, family=binomial)
summary(main_rare)

coef_main_rare <- (coeftest(main_rare, vcov=function(x)
  vcovHC(x, method= "arellano", cluster="group", type="HC1")))

# all variables regression (logit) column 5

main_glm <- glm(jurisdiction ~ nonwhite + rule_law +
                      sharia + capability + troops + 
                      nato + affinity + veto + cold_war, 
                    data = SOFA_models, family=binomial)
summary(main_glm)

coef_main_glm <- (coeftest(main_glm, vcov=function(x)
  vcovHC(x, method= "arellano", cluster="group", type="HC1")))


#################################### Appendix ########################################

# Table A4: Non-White Majority Host on U.S. SOFA Jurisdiction

# judicial independence instead of rule of law as stricter test

# all variables regression (rare logit) - column 6

main_jud_rare <- bayesglm(jurisdiction ~ nonwhite + judicial_ind +
                            sharia + capability + troops + 
                            nato + affinity + veto + cold_war, 
                          data = SOFA_models, family=binomial)
summary(main_jud_rare)

coef_main_jud_rare <- (coeftest(main_jud_rare, vcov=function(x)
  vcovHC(x, method= "arellano", cluster="group", type="HC1")))


# all variables regression (logit) - column 7

main_jud_glm <- glm(jurisdiction ~ nonwhite + judicial_ind +
                      sharia + capability + troops + 
                      nato + affinity + veto + cold_war, 
                    data = SOFA_models, family=binomial)
summary(main_jud_glm)

coef_main_jud_glm <- (coeftest(main_jud_glm, vcov=function(x)
  vcovHC(x, method= "arellano", cluster="group", type="HC1")))

# Table A5: Odds Ratios

exp_1 <- (exp(coef(coef_main_rare)))
exp_2 <- (exp(coef(coef_main_glm)))
exp_3 <- (exp(coef(coef_main_jud_rare)))
exp_4 <- (exp(coef(coef_main_jud_glm)))


# Correlation matrix 

correlation_data <- SOFA_models %>% 
  dplyr::select(jurisdiction, nonwhite, cold_war, sharia, nato,
                veto, rule_law, judicial_ind, capability, affinity, troops) 

core_table <- cor(correlation_data, use = "complete.obs") 


